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|  ABSTRACT 

A  method  of  subpixel  image  registration  is  proposed  that  employes  a  model  of  the 
correlation  surface  in  the  vicinity  of  the  registration  point. 
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1.  Introduction 


Registration  (or  matching  or  correlation)  is  the  process  of  determining  the 
geometric  transformation  which  is  required  in  order  to  conform  a  small  image 
with  a  large  image;  see  [lj.  The  simplest  transformation  is  pure  translation, 
where  the  matching  determines  the  translational  offset  of  the  small  image  in  the 
large  one.  For  each  possible  offset  one  defines  an  object  function,  whose 
extremum  (maximum  or  minimum  depending  on  the  function)  defines  the  desired 
matching.  For  this  area  (or  template)  correlation,  in  contrast  to  feature-space 
correlation,  the  object  function,  e.g.,  cross-correlation  or  sum  of  absolute 
differences,  directly  utilizes  the  pixel  (gray  level)  values  of  the  two  images.  The 
object  function  is  defined  only  for  a  discrete  grid  of  offsets  between  the  two 
images,  whereas  the  true  offset  of  the  registration  may  lie  somewhere  in  between 
these  grid  points.  If  the  geometric  distortion  between  the  two  images  is 
sufficiently  small,  subpixel  accuracy  in  the  registration  can  be  achieved. 

In  order  to  achieve  subpixel  registration  some  form  of  interpolation  must  be 
used.  A  straightforward  method  is  to  fit  some  analytic  surface  to  the  object 
function  samples  in  the  vicinity  of  the  extremum,  and  then  evaluate  analytically 
the  location  of  the  extremum  of  that  surface;  see  [2,  3].  A  similar  method, 
applied  by  the  Canadian  Center  for  Remote  Sensing,  is  to  compute  the  registra¬ 
tion  point  as  a  weighted  sum  of  neighboring  object  functions;  see  [4],  where  other 
methods,  such  as  Fourier  techniques,  are  also  mentioned. 

The  problem  with  these  methods  is  that  they  do  not  employ  a  model  of  the 
true  surface  of  the  object  function  in  the  neighborhood  of  the  registration  point. 


The  actual  surface  there  can  be  non-analytic  (cusped)  and  asymmetric  with 
respect  to  the  extremum  point.  In  what  follows  we  suggest  that  better  subpixel 
accuracy  in  the  registration  can  be  achieved  if  a  good  model  of  the  correlation 
surface  is  employed. 

2.  The  Method 

Our  main  idea,  then,  is  to  achieve  subpixel  accuracy  by  correlating  a  small 
grid  of  samples  of  the  object  function,  in  the  neighborhood  of  the  registration 
point,  with  a  model  of  the  correlation  surface  which  incorporates  its  (possible) 
non-analytic  and  asymmetric  nature.  The  steps  for  determining  the  subpixel 
registration  point  are  as  follows: 

(1.)  First,  the  approximate  (pixel)  registration  point,  P,  is  found  by  the 
usual  method:  The  object  function  between  the  small  image  and  the  large  image 
is  calculated  at  the  grid  points.  Then  the  point  P  where  this  function  achieves 
its  (global)  extremal  value  is  found. 

(2.)  Next,  a  model  for  the  correlation  surface  around  P  is  constructed  as 
follows:  A  patch  from  the  large  image  centered  at  P  and  of  the  size  of  the  small 
image  is  created.  Then  the  object  function  between  the  patch  and  the  large 
image  (the  “auto-object  function”)  is  calculated  on  a  small  grid  (say,  5  X  5  or 
7X7)  around  P .  The  model  surface  is  passed  through  these  ponts  and  is 
approximated  by  bilinear  or  spline  interpolation  between  them.  Note  that  the 
extremum  of  the  model  surface  is  precisely  at  P . 
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(3.)  Finally,  the  subpixel  registration  point  Q  is  found  by  correlating  a 
small  discrete  sample  (say  3  X  3  or  5  X  5  around  P )  of  the  object  function  cal¬ 
culated  in  (1),  with  the  continuous  model  surface  defined  in  (2).  The  mismatch 
between  the  grid  and  the  surface  defines  the  offset  Q  ,  the  subpixel  registration 
point,  and  P . 

If  there  exists  radiometric  (gray  level)  distortion  between  the  small  and  large 
images,  then  the  correlations  in  (1)  and  (3)  are  carried  out  after  normalization 
with  respect  to  the  mean  and  variance  values. 


3.  An  Example 

In  order  to  demonstrate  the  method  and  its  performance,  we  solve  below  a 
simple  one-dimensional  example.  The  “exact”  picture  is  assumed  to  be 

/(x)  =  sin(mi)  0  <  x  <  n  (1) 


where  m  is  some  integer.  The  large  image  consists  of  N  +  1  pixels  centered  at 


*  =  0,  1, - N 


with  gray  levels 


(2) 


ft  = 


sin(mx )  dx 
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cos  mh  (i  -  — )  -  cos  mh  (i  + 


(3) 
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The  small  image  is  created  by  sampling  the  following  continuous  function: 
g  (x  )  =  sin  (m  +  6m  )  (x  +  6x  )  +  e(x  ) 

(4) 

—  sin  (mi  +  m  6x  +  x  6m)  +  e(x  ) 

where  6m  introduces  a  geometric  distortion,  and  6x  introduces  a  subpixel  shift 
between  the  two  images.  e(x )  is  the  noise.  At  the  registration  point,  x  =  xreg, 
we  require  that,  without  noise,  /(xreg)  =  g(xreg),  i.e.,  xreg  6m  +  m  6x  =  0  and 


therefore 


<7  (x  )  =  sin  (m  +  6m  )  (x - —  x  )  +  e(x  ) 


The  small  image  consists  of  2  jw  +  1  pixels  centered  at 


where 


xj  ~  U  o  +  6j  +  j )  h  ,  ju 


^reg  {j  0  "b  ^3  )  ^ 


is  the  center  of  the  small  image.  Hence,  the  gray  levels  of  the  pixels  of  the  small 


image  are 


r  x,+i 

1]  =  J  g(x)  dx  =  — 
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Once  the  gray  levels  of  the  large  (/,  ),  and  small  [g }  )  images  are  given,  we 
can  calculate  the  cross-object  functions  defined  as  follows: 


h}'}]  [k  )  =  — —  Vj  \fk~j~9j 

”’<■  ;  = 


hfg‘]  (k  )  —  X]  fk-j  9j/  S  fk~j  S  9 j" 


j  = 


where  nu,  -  2 ,/„  +  1.  ju,  <  k  <  V  -  jw  . 

Next  we  create  the  auto-object  function  as  follows:  Assume  that  the  cross¬ 
object  function  defined  in  Eq.  (9)  assumes  its  global  extremal  value  (minimum  for 
and  maximum  for  hj^p)  at  the  grid  point  k  =  k*.  Then  the  auto-object 
functions  are  defined  as  follows: 


N'1’  (* )  =  E  I  A*,  -  /»•+  i  I 

"  J  = 

<*>-  E  A-,  A-,;/  E  Ah  E  As/ 

J  “Jw  J  “  Jw  J  Jw 


where  and  the  range  of  A;  are  defined  above,  and  the  limits  of  summation  are 
as  in  Eq.  (9).  Clearly  hj y  assumes  its  extremal  value  precisely  at  k  =  k*.  where 
h  {}  [k*]  =  0  and  h 1  ( k *)  —  1. 

The  subpixel  registration  can  now  be  found  by  the  two  methods  described  in 
Sections  1  and  2  above: 


Method  1 


Select  a  set  of  grid  points  (say  3  or  5)  around  k  —  k*  where  assumes  i t s 

global  extremal  value.  Fit  an  analytic  function  (say.  quadratic  or  quart ic 


"  *”*  ***  "'•*  •/  ^  ^ 


through  these  points.  Finally,  calculate  analytically  the  point  x  = 


w  he  re 


the  derivative  of  the  analytic  function  vanishes.  The  point  x  =  £re„  is  the  sub¬ 


pixel  registration  point,  and  the  fractional  error  due  to  the  tirst  met] 


sree  J  i 


>reg  *  reg 

h 


Method  2 


Here  the  registration  point  is  found  by  correlating  the  cross-object  function 


hf  ( k  ),  Eq,  (9),  with  the  auto-object  function  h,,  ( k  ),  Eq.  (10),  around  the  grid 


point  k  =  k‘.  We  first  continue  the  discrete  function  h,,  ( k  )  to  a  continuous 


function  of  s  whose  origin  is  located  at  k  =  k*  i.e.. 


(k  -  k*)  h 


where  h  is  the  pixel  size  defined  in  Eq.  (2).  The  continuous  function  h ,,  (.$ 


coincides  with  hjj  [k  )  at  he  grid  points,  and  is  interpolated  (say,  linearly) 


between  these  points. 


Next  we  construct  a  continuous  object  function  //  (s  )  between  hr,  (,->  )  and 


h  :rj  [k  )  defined  as  follows: 


//'"  Is 


=  —  E  I  hkl)(k*+J)  h"](*+jh) 


n>  j  =  -j> 


//l2)  (s) 


£  hj?Hk'+  j)h™(s  +  Jh)/D 
]  =  -  ik 


=  !  £  IV/'(e+  j))'  v  Hial(«  +  jh))'\ 


J  = 


i  = 


where  jh  —  1  or  2,  and  n;  =2  jh  +  1. 


The  registration  point  is  at  s  =  s'  where  H  (s  )  assumes  its  extremal  value 
(minimum  for  H{ ( s *)  and  maximum  for  H ^  ( s *)  ).  Since  H  [s  )  is  not  given  in 
closed  form,  s*  cannot  be  found  analytically,  but  has  to  be  evaluated  numerically 
(say,  by  Newton’s  method).  In  the  coordinate  system  of  the  large  image,  the  sub¬ 
pixel  registration  point  is  at  x  =  r/reg,  where 

/?reg  —  k  h  —  s  (14) 

and  the  fractional  error  due  to  the  second  method  is 

h -  (is) 
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